home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / bit / src / iprocess.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  23KB  |  905 lines
  1. /*
  2.  * $Id: iprocess.c,v 0.91 1994/02/20 00:53:18 zhao Pre-Release $
  3.  *
  4.  *. This file is part of BIT shareware package. After the two weeks of
  5.  *  free evaluation period, you are encouraged (required) to register
  6.  *  your copy for a small registration fee, which is $35 for personal use
  7.  *  and $50 for commercial, government and institutional use.
  8.  *
  9.  *  Copyright(c) 1993, 1994 by T.C. Zhao.
  10.  *  All rights reserved.
  11.  *
  12.  *  Permission to use, copy, and distribute this software in its entirety
  13.  *  for non-commercial purposes is hereby granted, provided that the
  14.  *  above shareware and copyright notices and this permission notice
  15.  *  appear in all copies and their documentation.
  16.  *
  17.  *  This software may be modified for your own use, but modified versions
  18.  *  may not be distributed without prior consent of the author.
  19.  *
  20.  *  This software is provided "as is" without expressed or implied
  21.  *  warranty of any kind.
  22.  *
  23.  *.
  24.  *
  25.  * Simple image processing routines. Not very isolated, but there
  26.  * are no-interactive tasks to be performed by any of these routines.
  27.  */
  28.  
  29. #if !defined(lint) && defined(F_ID)
  30. char *id_ip = "$Id: iprocess.c,v 0.91 1994/02/20 00:53:18 zhao Pre-Release $";
  31. #endif
  32.  
  33. #include "bit.h"
  34.  
  35. /********************************************************************
  36.  * General convlution routines.
  37.  * could be further optimized/hacked to go faster, but undoubtly
  38.  * will be uglier.
  39.  *******************************************************************/
  40.  
  41.  
  42. #define vectorP3(cm, m, c)                                           \
  43.         (cm[0]*m[c-1] + cm[1] * m[c] + cm[2]*m [c + 1])      \
  44.  
  45. #define vectorP5(cm, m, c)                                           \
  46.                (cm[0]*m[c-2]+cm[1]*m[c-1]+cm[2]*m[c]+                \
  47.                 cm[3]*m[c+1]+cm[4]*m[c+2])
  48.  
  49. /* Perform a 3x3 convolution at m(r,c)   */
  50.  
  51. #define conv3x3(cm,m,r,c)                                            \
  52.                (vectorP3(cm[0], m[r-1], c) +                         \
  53.                 vectorP3(cm[1], m[r], c) +                           \
  54.                 vectorP3(cm[2], m[r + 1], c))
  55.  
  56.  
  57. /* may require compiler flag -Wf,-XNh2000 to */
  58. #if 0
  59. #define conv5x5(cm,m,r,c)                                             \
  60.        (vectorP5(cm[0], m[r - 2], c) + vectorP5(cm[1], m[r - 1], c) + \
  61.         vectorP5(cm[2], m[r] , c) + vectorP5(cm[3], m[r + 1], c) +    \
  62.         vectorP5(cm[4], m[r + 2], c))
  63. #else
  64. static int
  65. conv5x5(register int **cm, register pc_t **m,
  66.     register int r, register int c)
  67. {
  68.     return (vectorP5(cm[0], m[r - 2], c) +
  69.         vectorP5(cm[1], m[r - 1], c) +
  70.         vectorP5(cm[2], m[r], c) +
  71.         vectorP5(cm[3], m[r + 1], c) +
  72.         vectorP5(cm[4], m[r + 2], c));
  73. }
  74. #endif
  75.  
  76. #define PC_check(pc)                                           \
  77.    do {                                                        \
  78.       if(pc < 0) pc = 0;                                       \
  79.       else if((pc /= weight) > PCMAXV)                         \
  80.           pc = PCMAXV;                                         \
  81.    } while (ZERO)
  82.  
  83. /*** thresholding ****/
  84.  
  85. #define TRSHD(new, old, t, i)                                 \
  86.      ((t && t[i])  ? ((Abs((int)(new)-(int)(old)) >= t[i]) ?  \
  87.           (new):(old)) : (new))
  88.  
  89. /****************************************************************
  90.  * Given RGB matrices and the convolution kernel, do the
  91.  * convolution.
  92.  ****************************************************************/
  93.  
  94. static int
  95. rgb_convolve(register pc_t **rm, register pc_t **gm, register pc_t **bm,
  96.          int nrow, int ncol, register int **mat, int cm_col, int cm_row,
  97.          register int *threshold, const char *what)
  98. {
  99.     register long weight = 0;
  100.     register int rnew, gnew, bnew;
  101.     register int row, col;
  102.     int i, j;
  103.     long rlines;
  104.  
  105.     /* get normalization factor */
  106.     for (row = 0; row < cm_row; row++)
  107.     for (col = 0; col < cm_col; col++)
  108.         weight += mat[row][col];
  109.  
  110.     if (weight <= 0)
  111.       {
  112.       Bark("Convolve", "Bad matrix weight: %ld", weight);
  113.       return -1;
  114.       }
  115.  
  116.     rlines = progress_report(what, nrow - cm_row / 2);
  117.  
  118.     if (cm_col == 3 && cm_row == 3)
  119.       {
  120.       nrow--;
  121.       ncol--;
  122.       for (row = 1; row < nrow; row++)
  123.         {
  124.         REPORT(row, rlines);
  125.         for (col = 1; col < ncol; col++)
  126.           {
  127.               rnew = conv3x3(mat, rm, row, col);
  128.               gnew = conv3x3(mat, gm, row, col);
  129.               bnew = conv3x3(mat, bm, row, col);
  130.               PC_check(rnew);
  131.               PC_check(gnew);
  132.               PC_check(bnew);
  133.               rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
  134.               gm[row][col] = TRSHD(gnew, gm[row][col], threshold, 1);
  135.               bm[row][col] = TRSHD(bnew, bm[row][col], threshold, 2);
  136.           }
  137.         }
  138.       }
  139.     else if (cm_col == 5 && cm_row == 5)
  140.       {
  141.       nrow -= 2;
  142.       ncol -= 2;
  143.       for (row = 2; row < nrow; row++)
  144.         {
  145.         REPORT(row, rlines);
  146.         for (col = 2; col < ncol; col++)
  147.           {
  148.               rnew = conv5x5(mat, rm, row, col);
  149.               gnew = conv5x5(mat, gm, row, col);
  150.               bnew = conv5x5(mat, bm, row, col);
  151.               PC_check(rnew);
  152.               PC_check(gnew);
  153.               PC_check(bnew);
  154.               rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
  155.               gm[row][col] = TRSHD(gnew, gm[row][col], threshold, 1);
  156.               bm[row][col] = TRSHD(bnew, bm[row][col], threshold, 2);
  157.           }
  158.         }
  159.       }
  160.     else
  161.       {
  162.       int ccol, ii, jj;
  163.       int cm_2_row = cm_row / 2;
  164.       int cm_2_col = cm_col / 2;
  165.       for (row = cm_2_row; row < nrow - cm_2_row; row++)
  166.         {
  167.         REPORT(row, rlines);
  168.         for (col = cm_2_col; col < ncol - cm_2_col; col++)
  169.           {
  170.               rnew = gnew = bnew = 0;
  171.               ccol = col - cm_2_col;
  172.               for (i = 0; i < cm_row; i++)
  173.             {
  174.                 ii = row - cm_2_row + i;
  175.                 for (jj = ccol, j = 0; j < cm_col; jj++, j++)
  176.                   {
  177.                   rnew += mat[i][j] * rm[ii][jj];
  178.                   gnew += mat[i][j] * gm[ii][jj];
  179.                   bnew += mat[i][j] * bm[ii][jj];
  180.                   }
  181.             }
  182.  
  183.               /* clamp to minmax */
  184.               PC_check(rnew);
  185.               PC_check(gnew);
  186.               PC_check(bnew);
  187.               rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
  188.               gm[row][col] = TRSHD(gnew, gm[row][col], threshold, 1);
  189.               bm[row][col] = TRSHD(bnew, bm[row][col], threshold, 2);
  190.           }
  191.         }
  192.       }
  193.     remove_progress_report();
  194.     return 0;
  195. }
  196.  
  197. /*****************************************************************
  198.  * Similar  to rgb_convolve except that the input is in grayscale.
  199.  * Note gm & bm are un-used
  200.  ******************************************************************/
  201.  
  202. /* ARGSUSED */
  203. static int
  204. gray_convolve(register pc_t **rm, pc_t **gm, pc_t **bm,
  205.           int nrow, int ncol, register int **mat, int cm_col, int cm_row,
  206.           register int *threshold, const char *what)
  207. {
  208.     register long weight = 0;
  209.     register int rnew;
  210.     register int row, col, i, j, ii, jj;
  211.     long rlines;
  212.  
  213.     /* get normalization factor */
  214.     for (row = 0; row < cm_row; row++)
  215.     for (col = 0; col < cm_col; col++)
  216.         weight += mat[row][col];
  217.  
  218.     if (weight <= 0)
  219.       {
  220.       Bark("Convolve", "Bad matrix weight: %ld", weight);
  221.       return -1;
  222.       }
  223.  
  224.     rlines = progress_report(what, nrow - cm_row / 2);
  225.     if (cm_col == 3 && cm_row == 3)
  226.       {
  227.       nrow--;
  228.       ncol--;
  229.       for (row = 1; row < nrow; row++)
  230.         {
  231.         REPORT(row, rlines);
  232.         for (col = 1; col < ncol; col++)
  233.           {
  234.               rnew = conv3x3(mat, rm, row, col);
  235.               PC_check(rnew);
  236.               rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
  237.           }
  238.         }
  239.       }
  240.     else if (cm_col == 5 && cm_row == 5)
  241.       {
  242.       nrow -= 2;
  243.       ncol -= 2;
  244.       for (row = 2; row < nrow; row++)
  245.         {
  246.         REPORT(row, rlines);
  247.         for (col = 2; col < ncol; col++)
  248.           {
  249.               rnew = conv5x5(mat, rm, row, col);
  250.               PC_check(rnew);
  251.               rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
  252.           }
  253.         }
  254.       }
  255.     else
  256.       {
  257.       int ccol;
  258.       int cm_2_row = cm_row / 2;
  259.       int cm_2_col = cm_col / 2;
  260.       for (row = cm_2_row; row < nrow - cm_2_row; row++)
  261.         {
  262.         REPORT(row, rlines);
  263.         for (col = cm_2_col; col < ncol - cm_2_col; col++)
  264.           {
  265.               rnew = 0;
  266.               ccol = col - cm_2_col;
  267.               for (i = 0; i < cm_row; i++)
  268.             {
  269.                 ii = row - cm_2_row + i;
  270.                 for (jj = ccol, j = 0; j < cm_col; jj++, j++)
  271.                   {
  272.                   rnew += mat[i][j] * rm[ii][jj];
  273.                   }
  274.             }
  275.  
  276.               /* clamp to minmax */
  277.               PC_check(rnew);
  278.               rm[row][col] = TRSHD(rnew, rm[row][col], threshold, 0);
  279.           }
  280.         }
  281.       }
  282.     remove_progress_report();
  283.     return 0;
  284. }
  285.  
  286. /*************************************************************
  287.  * Global interface to convolution routines
  288.  *
  289.  * do the convolution on an image or a part of it.
  290.  *************************************************************/
  291.  
  292. int
  293. img_convolv(IPTR im, register int **m, int nrow, int ncol,
  294.         int *t, const Rect_t * sr, const char *msg)
  295. {
  296.     register pc_t **rm, **gm, **bm;
  297.     int status;
  298.     int total = sr->h * sr->w;
  299.     rgba_t **subi;
  300.  
  301.     /* always demand proper region */
  302.     if (sr->w < ncol || sr->h < nrow)
  303.     return -1;
  304.  
  305.     fl_deactivate_all_forms();
  306.  
  307.     /*
  308.      * smoothing can only be done to RGB images. Turn maped gray image to
  309.      * gray and all others to RGB
  310.      */
  311.  
  312.     if (!IS_CPACK(im))
  313.       {
  314.       (void) img_convert_type(im, IS_GRAY(im) ? T_GRAY : T_RGBA);
  315.  
  316.       /* re-display */
  317.       im->io->display(im, -1, 0);
  318.       }
  319.  
  320.     /* check if to be done on the entire image */
  321.  
  322.     if (!equal_rect(sr, img_rect(im)))
  323.       {
  324.  
  325.       /* not whole image */
  326.       subi = get_subimage(im, sr->x, sr->y, sr->w, sr->h);
  327.  
  328.       /* gm & bm are un-used n convolution is grayscale image */
  329.       rm = get_mat(sr->h, sr->w, sizeof(pc_t));
  330.       gm = get_mat(sr->h, sr->w, sizeof(pc_t));
  331.       bm = get_mat(sr->h, sr->w, sizeof(pc_t));
  332.  
  333.       /* unpack RGBA to R, G, B. */
  334.       unpack_mat(subi[0], rm[0], gm[0], bm[0], total);
  335.       }
  336.     else
  337.       {
  338.       subi = im->mraster;
  339.       img_get_rgb(im);
  340.       rm = get_RM(im);
  341.       gm = get_GM(im);
  342.       bm = get_BM(im);
  343.       }
  344.  
  345.  
  346.     status = (IS_GRAY(im) ? gray_convolve : rgb_convolve)
  347.     (rm, gm, bm, sr->h, sr->w, m, nrow, ncol, t, msg);
  348.  
  349.     /* pack_mat detects gray & rgba by checking if R == G */
  350.     pack_mat(subi[0], rm[0], IS_GRAY(im) ? rm[0] : gm[0], bm[0], total);
  351.  
  352.     /* change the image in core */
  353.     put_subimage(im, subi, sr, 1);
  354.  
  355.     /* change the image in framebuffer */
  356.     set_current_window(win_id);
  357.     Rectwrite(sr->x, sr->y, sr->x + sr->w - 1, sr->y + sr->h - 1, subi[0]);
  358.  
  359.     /* do the clean up */
  360.     if (subi == im->mraster)
  361.       {
  362.       img_free_rgbmem(im);
  363.       }
  364.     else
  365.       {
  366.       free_mat(rm);
  367.       free_mat(gm);
  368.       free_mat(bm);
  369.       free_mat(subi);
  370.       }
  371.     fl_activate_all_forms();
  372.  
  373.     return status;
  374. }
  375.  
  376.  
  377. /*
  378.  * sum of two derivatives and convert to colormapped image only looks ok
  379.  * followed by a thresholding with coledit
  380.  */
  381. #define get_derivative(r0, r1, r2, i)                            \
  382.    do {                                                          \
  383.       d1 = (int) r0[i + 1] - (int) r0[i - 1] +                   \
  384.            2 * ((int) r1[i + 1] - (int) r1[i - 1]) +             \
  385.            (int) r2[i + 1] - (int) r2[i - 1];                    \
  386.       d2 = (int) r2[i - 1] + 2 * (int) r2[i] + (int) r2[i + 1] - \
  387.            ((int) r0[i - 1] + 2 * (int) r0[i] + r0[i + 1]);      \
  388.    } while (ZERO)
  389.  
  390. int
  391. img_edge(IPTR im)
  392. {
  393.     int i, j;
  394.     register pc_t *r0, *r1, *r2;
  395.     register pc_t *g0, *g1, *g2;
  396.     register pc_t *b0, *b1, *b2;
  397.     long rline;
  398.     ci_t **m;
  399.     int d1, d2, d;
  400.     int gray = IS_GRAY(im);
  401.  
  402.  
  403.     show_busy("Just a moment...");
  404.     if (img_get_rgb(im) < 0)
  405.     return -1;
  406.  
  407.     if (!IS_CI(im))
  408.       {
  409.       /* free the old image raster */
  410.       img_free_rastermem(im);
  411.  
  412.       /* new image */
  413.       m = get_mat(im->h, im->w, sizeof(ci_t));
  414.       }
  415.     else
  416.       {
  417.       m = im->mraster;
  418.       }
  419.  
  420.     rline = progress_report("Searching ...", im->h);
  421.     memset(m[0], 0, sizeof(ci_t) * im->w);
  422.     memset(m[im->h - 1], 0, sizeof(ci_t) * im->w);
  423.  
  424.     for (j = 0; j < im->h; j++)
  425.     m[j][0] = m[j][im->w - 1] = 0;
  426.  
  427.     im->cmap->colors = PCMAXV;    /* final range */
  428.  
  429.     if (gray)
  430.       {
  431.       for (j = 1; j < im->h - 1; j++)
  432.         {
  433.         r0 = get_RM(im)[j - 1];
  434.         r1 = get_RM(im)[j];
  435.         r2 = get_RM(im)[j + 1];
  436.         REPORT(j, rline);
  437.         for (i = 1; i < im->w - 1; i++)
  438.           {
  439.               get_derivative(r0, r1, r2, i);
  440.               d = ((float) (Abs(d1) + Abs(d2)) / 2.6);
  441.               if (d > im->cmap->colors)
  442.               d = im->cmap->colors;
  443.               m[j][i] = d;
  444.           }
  445.         }
  446.       }
  447.     else
  448.       {
  449.       for (j = 1; j < im->h - 1; j++)
  450.         {
  451.         r0 = get_RM(im)[j - 1];
  452.         r1 = get_RM(im)[j];
  453.         r2 = get_RM(im)[j + 1];
  454.         g0 = get_GM(im)[j - 1];
  455.         g1 = get_GM(im)[j];
  456.         g2 = get_GM(im)[j + 1];
  457.         b0 = get_BM(im)[j - 1];
  458.         b1 = get_BM(im)[j];
  459.         b2 = get_BM(im)[j + 1];
  460.         REPORT(j, rline);
  461.         for (i = 1; i < im->w - 1; i++)
  462.           {
  463.               get_derivative(r0, r1, r2, i);
  464.               d = Abs(d1) + Abs(d2);
  465.               get_derivative(g0, g1, g2, i);
  466.               d += Abs(d1) + Abs(d2);
  467.               get_derivative(b0, b1, b2, i);
  468.               d += Abs(d1) + Abs(d2);
  469.               d /= 8;
  470.               if (d > im->cmap->colors)
  471.               d = im->cmap->colors;
  472.               m[j][i] = d;
  473.           }
  474.         }
  475.       }
  476.     img_free_rgbmem(im);
  477.     /*
  478.      * probably should get the map that matches the original image closest in
  479.      * terms warmth
  480.      */
  481.     get_bit_defmap(!gray, im->cmap->ct[0],
  482.            im->cmap->ct[1],
  483.            im->cmap->ct[2],
  484.            im->cmap->colors);
  485.     remove_progress_report();
  486.     return fill_image_struct(im, m, im->h, im->w, gray ? T_GMAP : T_CMAP);
  487. }
  488.  
  489.  
  490. /********************************************************************
  491.  * histogram equalization with thresholding(not written yet)
  492.  ********************************************************************/
  493.  
  494. #include "lookup.h"
  495.  
  496. int
  497. img_enhance(IPTR im)
  498. {
  499.     register pc_t *r, *g, *b;
  500.     register int row, col, i;
  501.     double fact;
  502.     long rlines, total = im->h * im->w;
  503.     register unsigned int *hist = im->cmap->p_h.hist[3];
  504.     static long sum[PCMAX];
  505.  
  506.     if (img_get_rgb(im) < 0)
  507.     return -1;
  508.  
  509.     img_get_histogram(im);
  510.     memset(sum, 0, sizeof(sum));
  511.     sum[0] = hist[0];
  512.  
  513.     for (i = 1; i < PCMAX; i++)
  514.     sum[i] = sum[i - 1] + hist[i];
  515.  
  516.     /* get the mapping function */
  517.  
  518.     fact = (double) (PCMAX - 1) / total;
  519.     for (i = 0; i < PCMAX; i++)
  520.     sum[i] = (fact * sum[i]);
  521.  
  522.     rlines = progress_report("Enhancing Contrast ...", im->h);
  523.     r = get_RM(im)[0];
  524.     g = get_GM(im)[0];
  525.     b = get_BM(im)[0];
  526.  
  527.     for (row = 0; row < im->h; row++)
  528.       {
  529.       REPORT(row, rlines);
  530.       for (col = 0; col < im->w; col++, r++, g++, b++)
  531.         {
  532.         *r = sum[*r];
  533.         *g = sum[*g];
  534.         *b = sum[*b];
  535.         }
  536.       }
  537.     remove_progress_report();
  538.     row = img_rgb_to_cpack(im);
  539.     img_free_rgbmem(im);
  540.     return row;
  541. }
  542.  
  543. /*********************************************************************
  544.  * normalize a image. Similar to enhance, the normal curve will based
  545.  * on gray levels and comman to all.
  546.  **********************************************************************/
  547. int
  548. img_norm(IPTR im)
  549. {
  550.     register pc_t *r, *g, *b, *rs;
  551.     pc_t lut[PCMAX + 1];
  552.     static int crit = 10, hi;
  553.     int i;
  554.  
  555.     if (getint("Enter threshold", &crit, 0, (PCMAX - 1) / 2, 0) < 0)
  556.     return -1;
  557.  
  558.     if (img_get_rgb(im) < 0)
  559.     return -1;
  560.  
  561.     hi = PCMAX - 1 - crit;
  562.     for (i = 0; i < crit; i++)
  563.     lut[i] = 0;
  564.     for (i = hi; i < PCMAX; i++)
  565.     lut[i] = PCMAX - 1;
  566.     for (i = crit; i < hi; i++)
  567.     lut[i] = (hi * (i - crit) / (float) (hi - crit));
  568.  
  569.     r = get_RM(im)[0];
  570.     g = get_GM(im)[0];
  571.     b = get_BM(im)[0];
  572.  
  573.     for (rs = r + im->w * im->h; r < rs; r++, g++, b++)
  574.       {
  575.       *r = lut[*r];
  576.       *g = lut[*g];
  577.       *b = lut[*b];
  578.       }
  579.  
  580.     hi = img_rgb_to_cpack(im);
  581.     img_free_rgbmem(im);
  582.     return hi;
  583. }
  584.  
  585. /***************************************************************
  586.  * Show (RGB seperated) histograms of the current image.
  587.  * Note the histogram has two extra entries for max and ave
  588.  ***************************************************************/
  589. int
  590. img_show_histogram(IPTR im)
  591. {
  592.     int i, j;
  593.     Hist_t *hist;
  594.     float x[PCMAX], y[4][PCMAX], *yy[4], totalpix;
  595.     char ave[4][50];
  596.  
  597.     if (!image_ready(im, "Histogram"))
  598.     return -1;
  599.  
  600.     if (img_get_histogram(im) < 0)
  601.     return -1;
  602.  
  603.     set_chart4_help(HELP_HIST);
  604.  
  605.     hist = im->cmap->p_h.hist;
  606.     totalpix = im->h * im->w;
  607.  
  608.     for (i = 0; i < PCMAX; i++)
  609.       {
  610.       x[i] = i;
  611.       for (j = 0; j < 4; j++)
  612.           y[j][i] = (double) hist[j][i] / totalpix;
  613.       }
  614.  
  615.     for (i = 0; i < 4; i++)
  616.       {
  617.       yy[i] = y[i];
  618.       sprintf(ave[i], "Ave=%u Peak=%.1g at %u",
  619.       (unsigned) hist[i][PCMAX + 1], (double) hist[i][PCMAX] / totalpix,
  620.           hist[i][PCMAX + 2]);
  621.       }
  622.     set_chart4_text(ave[0], ave[1], ave[2], ave[3]);
  623.     show_chart4("Histogram", x, yy, PCMAX);
  624.     return 0;
  625. }
  626.  
  627. /* ARGSUSED */
  628. int
  629. img_fft(IPTR im)
  630. {
  631.     to_be_written("FFT");
  632.     return 0;
  633. }
  634.  
  635. /******************************************************************
  636.  * Compute pixel-pixel difference of two identically sized images.
  637.  * And show the resulting difference as an image
  638.  *******************************************************************/
  639. int
  640. img_diff(IPTR im1)
  641. {
  642.     IPTR im2;
  643.     register rgba_t *ras2, *ras, *rs;
  644.     register int r1, r2, r, g1, g2, g, b1, b2, b;
  645.     const char *fname = getfilename("Reference Image", ".", "*.*", "", 1);
  646.  
  647.     if (!fname || !*fname)
  648.     return -1;
  649.  
  650.     if ((im2 = open_image(fname)))
  651.       {
  652.       close_image(im2);
  653.  
  654.       /* currently JPEG desc fakes image size to 1 */
  655.       if ((im2->w != 1) && (im2->w != im1->w) && im1->h != im2->h)
  656.         {
  657.         Bark("IDiff", "%s & %s are of different size", fname,
  658.              im1->ifile);
  659.         free_image(im2);
  660.         return -1;
  661.         }
  662.  
  663.       free_image(im2);
  664.       im2 = load_image(fname);
  665.       }
  666.  
  667.     if (!im2 || !im2->ok)
  668.     return -1;
  669.  
  670.     if (im1->w != im2->w || im1->h != im2->h)
  671.       {
  672.       Bark("IDiff", "%s and %s are of different size", fname, im1->ifile);
  673.       free_image(im2);
  674.       return -1;
  675.       }
  676.  
  677.     strcpy(im1->ifile, "Idiff");
  678.  
  679.     /* both image should be in RGB */
  680.     if (img_convert_type(im1, T_RGBA) < 0 ||
  681.     img_convert_type(im2, T_RGBA))
  682.       {
  683.       Bark("IDiff", "type conversion failed");
  684.       free_image(im2);
  685.       return -1;
  686.       }
  687.  
  688.     /* get the difference: put the results into IM1 */
  689.     ras = ((rgba_t **) im1->mraster)[0];
  690.     ras2 = ((rgba_t **) im2->mraster)[0];
  691.     for (rs = ras + im1->w * im1->h; ras < rs; ras++, ras2++)
  692.       {
  693.       Unpack(*ras, r1, g1, b1);
  694.       Unpack(*ras2, r2, g2, b2);
  695.       r = Abs(r1 - r2);
  696.       g = Abs(g1 - g2);
  697.       b = Abs(b1 - b2);
  698.       *ras = Pack(r, g, b);
  699.       }
  700.     free_image(im2);
  701.     return 0;
  702. }
  703.  
  704. /******************************************************************
  705.  * Median filtering on a 3x3 grid
  706.  ***************************************************************{**/
  707.  
  708. /****************************************************************
  709.  * Given RGB matrices and the convolution kernel, do the
  710.  * convolution.
  711.  ****************************************************************/
  712.  
  713. static int med3x3(pc_t *, pc_t *, pc_t *);
  714.  
  715. static int
  716. rgb_median(register pc_t **rm, register pc_t **gm, register pc_t **bm,
  717.        int nrow, int ncol, int cm_row, int cm_col)
  718. {
  719.     register int row, col;
  720.     long rlines;
  721.  
  722.  
  723.     rlines = progress_report("Med3x3", nrow - cm_row / 2);
  724.  
  725.     if (cm_col == 3 && cm_row == 3)
  726.       {
  727.       nrow--;
  728.       ncol--;
  729.       for (row = 1; row < nrow; row++)
  730.         {
  731.         REPORT(row, rlines);
  732.         for (col = 1; col < ncol; col++)
  733.           {
  734.               rm[row][col] = med3x3(rm[row - 1] + col - 1,
  735.                         rm[row] + col - 1,
  736.                         rm[row + 1] + col - 1);
  737.               gm[row][col] = med3x3(gm[row - 1] + col - 1,
  738.                         gm[row] + col - 1,
  739.                         gm[row + 1] + col - 1);
  740.               bm[row][col] = med3x3(bm[row - 1] + col - 1,
  741.                         bm[row] + col - 1,
  742.                         bm[row + 1] + col - 1);
  743.           }
  744.         }
  745.       }
  746.     remove_progress_report();
  747.     return 0;
  748. }
  749.  
  750. static int
  751. gray_median(register pc_t **rm, register pc_t **gm, register pc_t **bm,
  752.         int nrow, int ncol, int cm_row, int cm_col)
  753. {
  754.     register int row, col;
  755.     long rlines;
  756.  
  757.  
  758.     rlines = progress_report("Med3x3", nrow - cm_row / 2);
  759.  
  760.     if (cm_col == 3 && cm_row == 3)
  761.       {
  762.       nrow--;
  763.       ncol--;
  764.       for (row = 1; row < nrow; row++)
  765.         {
  766.         REPORT(row, rlines);
  767.         for (col = 1; col < ncol; col++)
  768.           {
  769.               rm[row][col] = med3x3(rm[row - 1] + col - 1,
  770.                         rm[row] + col - 1,
  771.                         rm[row + 1] + col - 1);
  772.           }
  773.         }
  774.       }
  775.     remove_progress_report();
  776.     return 0;
  777. }
  778.  
  779. /**********************************************************************
  780.  * Perform median filtering on the part of the image, bounded by mr
  781.  **********************************************************************/
  782.  
  783. int
  784. img_median_filter(IPTR im, Rect_t * mr)
  785. {
  786.     rgba_t **subi;        /* subimage */
  787.     pc_t **rm, **bm, **gm;
  788.     int status;
  789.     int total = mr->w * mr->h;
  790.  
  791.     /* must be done to a region larger  3x3 */
  792.     if (mr->w < 3 || mr->h < 3)
  793.     return -1;
  794.  
  795.     /* we can only do this on RGB images */
  796.     if (!IS_CPACK(im))
  797.       {
  798.       img_convert_type(im, IS_GRAY(im) ? T_GRAY : T_RGBA);
  799.  
  800.       /* reset display mode and show image */
  801.       im->io->display(im, -1, 0);
  802.       }
  803.  
  804.     /* check if on the entire image */
  805.     if (!equal_rect(mr, img_rect(im)))
  806.       {
  807.  
  808.       if (!(subi = get_subimage(im, mr->x, mr->y, mr->w, mr->h)))
  809.           return -1;
  810.  
  811.       rm = get_mat(mr->h, mr->w, sizeof(**rm));
  812.       gm = get_mat(mr->h, mr->w, sizeof(**gm));
  813.       bm = get_mat(mr->h, mr->w, sizeof(**bm));
  814.  
  815.       /* unpack RGBA to R, G, B. */
  816.       unpack_mat(subi[0], rm[0], gm[0], bm[0], total);
  817.       }
  818.     else
  819.       {
  820.       img_get_rgb(im);
  821.       subi = im->mraster;
  822.       rm = get_RM(im);
  823.       gm = get_GM(im);
  824.       bm = get_BM(im);
  825.       }
  826.  
  827.     /* now do the stuff */
  828.     status = (IS_GRAY(im) ? gray_median : rgb_median)
  829.     (rm, gm, bm, mr->h, mr->w, 3, 3);
  830.  
  831.     /* pack_mat detects gray & rgba by checking if R == G */
  832.     pack_mat(subi[0], rm[0], IS_GRAY(im) ? rm[0] : gm[0], bm[0], total);
  833.  
  834.     /* change the image in core */
  835.     put_subimage(im, subi, mr, 1);
  836.  
  837.     /* change the image in framebuffer */
  838.     set_current_window(win_id);
  839.     Rectwrite(mr->x, mr->y, mr->x + mr->w - 1, mr->y + mr->h - 1, subi[0]);
  840.  
  841.     /* clean up */
  842.     if (rm != get_RM(im))
  843.       {
  844.       free_mat(rm);
  845.       free_mat(gm);
  846.       free_mat(bm);
  847.       free_mat(subi);
  848.       }
  849.     else
  850.       {
  851.       img_free_rgbmem(im);
  852.       }
  853.  
  854.     return status;
  855. }
  856.  
  857.  
  858. /*
  859.  * Median Finding on a 3-by-3 Grid
  860.  * by Alan Paeth
  861.  * from "Graphics Gems", Academic Press, 1990
  862. */
  863.  
  864. #define s2(a,b) if ((t = b-a) < 0) {a += t; b -= t;}
  865. #define mn3(a,b,c) s2(a,b); s2(a,c)
  866. #define mx3(a,b,c) s2(b,c); s2(a,c)
  867. #define mnmx3(a,b,c) mx3(a,b,c); s2(a,b)
  868. #define mnmx4(a,b,c,d) s2(a,b); s2(c,d); s2(a,c); s2(b,d)
  869. #define mnmx5(a,b,c,d,e) s2(a,b); s2(c,d); mn3(a,c,e); mx3(b,d,e)
  870. #define mnmx6(a,b,c,d,e,f) s2(a,d); s2(b,e); s2(c,f);\
  871.                             mn3(a,b,c); mx3(d,e,f)
  872. /*
  873.  * Find median on a 3x3 input box of integers.
  874.  * b1, b2, b3 are pointers to the left-hand edge of three
  875.  * parallel scan-lines to form a 3x3 spatial median.
  876.  * Rewriting b2 and b3 as b1 yields code which forms median
  877.  * on input presented as a linear array of nine elements.
  878.  */
  879. static int
  880. med3x3(pc_t *b1, pc_t *b2, pc_t *b3)
  881. {
  882.     register int r1, r2, r3, r4, r5, r6;
  883.     register int t;
  884.  
  885.     r1 = *b1++;
  886.     r2 = *b1++;
  887.     r3 = *b1++;
  888.     r4 = *b2++;
  889.     r5 = *b2++;
  890.     r6 = *b2++;
  891.  
  892.     mnmx6(r1, r2, r3, r4, r5, r6);
  893.     r1 = *b3++;
  894.     mnmx5(r1, r2, r3, r4, r5);
  895.     r1 = *b3++;
  896.     mnmx4(r1, r2, r3, r4);
  897.     r1 = *b3++;
  898.     mnmx3(r1, r2, r3);
  899.     return (r2);
  900. }
  901.  
  902. /***********************************************************************
  903.  * End of Median filter
  904.  *********************************************************************}*/
  905.